///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file linear_method.cpp
///
/// \brief   implementation of the fast linear method wavefunction optimizer
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#include<src/pcps.h>
#include<src/hamiltonian.h>
#include<src/wavefunction.h>
#include<src/vmc_walker.h>
#include<src/accumulator.h>
#include<src/configuration.h>
#include<src/davidson.h>
#include<src/deriv_ratio_accumulator.h>

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   optimizes the wavefunction using the fast linear method
///
/// \param[in]      userinp  the user's input options
/// \param[in]      ham      the Hamiltonian
/// \param[in]      wfn      the wavefunction
///
/// \return the compound index
///
///////////////////////////////////////////////////////////////////////////////////////////////////
template <class T> void pcps::linear_method_optimization(const pcps::Input & userinp,
                                                         const pcps::Hamiltonian<T> & ham,
                                                         pcps::Wavefunction<T> & wfn) {

  // just temporary to force the configuration class to instantiate
  pcps::Configuration<T> config(userinp);

//  // get MPI info
//  const int nproc  = pcps::mpi::size();
//  const int myrank = pcps::mpi::rank();
//
//  // iteratively update the wavefunction
//  for (int iter = 0; iter < userinp.max_macro_loop(); iter++) {
//
//    // prepare a walker for use in sampling
//    pcps::VMC_Walker<T> walker(userinp);
//
//    // set the guiding function as the overall wavefunction coefficient
//    if ( wfn.overall_coefficient() )
//      walker.gf() = wfn.overall_coefficient();
//    else
//      throw pcps::Exception("wavefunction does not have an overall coefficient object in pcps::linear_method_optimization");
//
//    // burn in the markov chain
//    pcps::sample(userinp, userinp.burn_in_length(), walker);
//
//    // save the configuration after burn in
//    const pcps::Configuration<T> saved_config(walker.config());
//
//    // save the state of the random number generator after burn in
//    const pcps::LaggedFibonacci<> saved_lfg(walker.lfg());
//
//    // get the number of wavefunction variables
//    const int n = 1 + wfn.overall_coefficient()->nfvar();
//
//    // create an accumulator for value and derivative ratios
//    boost::shared_ptr<pcps::Accumulator<T> > dra( new pcps::DerivRatioAccumulator<T>(userinp,
//                                                                                     wfn.overall_coefficient(),
//                                                                                     wfn.overall_coefficient()) );
//
//    // accumulate the value and derivative ratios
//    walker.accumulators().push_back(dra);
//    pcps::sample(userinp, userinp.sample_length_per_proc(), walker);
//    walker.accumulators().clear();
//
//    // create the object that will multiply by the Hamiltonian
//    pcps::HamiltonianActor<T> h_actor(userinp, ham, wfn, n, saved_config, saved_lfg);
//
//    // create the object that will multiply by the overlap matrix
//    pcps::OverlapActor<T> s_actor(userinp, n, saved_config, saved_lfg, dra);
//
//    // declare variable to hold the energy from the davidson iterations
//    T E;
//
//    // create the object that will multiply by the preconditioner
//    pcps::LinearMethodPreconActor<T> precon(userinp, n, saved_config, saved_lfg, dra, E);
//
//    // allocate a vector to hold the davidson solution and put the initial guess (the current wavefunction) in it
//    std::vector<T> x(n, pcps::zero<T>());
//    x.at(0) = pcps::unity<T>();
//
//    // use the davidson solver to find the new wavefunction
//    pcps::davidson(n, userinp.max_krylov_iter(), userinp.krylov_thresh(), E, &x.at(0), h_actor, s_actor, precon);
//
//    // get the wavefunction variables in vector form
//    std::vector<T> y(n-1, pcps::zero<T>());
//    wfn.to_vector(y);
//
//    // add the update
//    pcps::xaxpy(n, 1.0 / x.at(0), &x.at(1), 1, &y.at(0), 1);
//
//    // broadcast the wavefunction variables to all processes
//    pcps::mpi::bcast(y);
//
//    // place the updated variables in the wavefunction
//    wfn.from_vector(y);
//
//  }

}


template void pcps::linear_method_optimization(const pcps::Input &,
                                               const pcps::Hamiltonian<double> &,
                                               pcps::Wavefunction<double> &);
template void pcps::linear_method_optimization(const pcps::Input &,
                                               const pcps::Hamiltonian<std::complex<double> > &,
                                               pcps::Wavefunction<std::complex<double> > &);
